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ABSTRACT 

Aims. The transiting planet WASP-12 b was identified as a potential target for transit timing studies because a departure from a linear ephemeris 
was reported in the literature. Such deviations could be caused by an additional planet in the system. We attempt to confirm the existence of 
claimed variations in transit timing and interpret its origin. 

Methods. We organised a multi-site campaign to observe transits by WASP-12 b in three observing seasons, using 0.5-2.6-metre telescopes. 
Results. We obtained 61 transit light curves, many of them with sub-millimagnitude precision. The simultaneous analysis of the best-quality 
datasets allowed us to obtain refined system parameters, which agree with values reported in previous studies. The residuals versus a linear 
ephemeris reveal a possible periodic signal that may be approximated by a sinusoid with an amplitude of 0.00068 ± 0.00013 d and period of 
500 ± 20 orbital periods of WASP-12 b. The joint analysis of timing data and published radial velocity measurements results in a two-planet model 
which better explains observations than single-planet scenarios. We hypothesize that WASP-12 b might be not the only planet in the system and 
there might be the additional 0. 1 Mj up body on a 3.6-d eccentric orbit. A dynamical analysis indicates that the proposed two-planet system is stable 
over long timescales. 
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1. Introduction 

The transiting extrasolar planet WASP-12 b was found to be one 
of the most intensely irradiated planets (Hebb et al. 12009b . It 
orbits a G-type star on a tight orbit with a semimajor axis of 
a h = 0.0229 ± 0.0008 AU and an orbital period of only -1.09 d. 
The planet's proximity to the star results in a high equilibrium 
temperature of 2500 K, thus inducing numerous studies of the 
properties of the planetary atmosphere (Li et al. 120101 Lai et 
al. 120101 Fossati et al. |2010bl Croll et al. I20TT1 Crossfield et 
al. 120121) . In the discovery paper, the planet was found to have 
a mass of m\, = 1.41 + 0.10 Mj up and an unexpectedly large 
radius of R\, = 1.79 ± 0.09 Rj up . The host star was found to 
have an effective temperature of T e g = 6300^^ K and a mass 
of M, = 1.35 + 0.14 M Q . Follow-up observations, primarily 
high-precision photometry, confirmed the large planetary radius 
(Maciejewski et al. 1201 II Chan et al. 1201 Q . The system param- 



* Partly based on (1) data collected with the Nordic Optical 
Telescope, operated on the island of La Palma jointly by Denmark, 
Finland, Iceland, Norway, and Sweden, in the Spanish Observatorio del 
Roque de los Muchachos of the Instituto de Astrofisica de Canarias, 
(2) observations made at the Centro Astronomico Hispano Aleman 
(CAHA), operated jointly by the Max-Planck Institut fur Astronomie 
and the Instituto de Astroffsica de Andalucia (CSIC), (3) data collected 
with telescopes at the Rozhen National Astronomical Observatory, and 
(4) observations obtained with telescopes of the University Observatory 
lena, which is operated by the Astrophysical Institute of the Friedrich- 
Schiller-University. 



eters were also refined by Southworth (120121 ) who used high- 
precision literature data in a joint analysis. 

The orbital eccentricity of WASP-12 b, e b = 0.049 ± 0.015, 
was initially determined from the radial velocity (RV) mea- 
surements (Hebb et al. 2009). Lopez-Morales et al. (2010) ob- 
served occultations whose timing supports a non-circular or- 
bit. On the other hand, further observations of occultations by 
Campo et al. (2011) yielded mid-occultation times at orbital 
phases 0.5010 + 0.0006 and 0.5006 + 0.0007, consistent with 
a circular orbit. The same conclusion was reached by Husnoo 
et al. (1201 11 ) who analysed spectroscopic data enhanced by new 
RV measurements. Croll et al. ( 1201 II ) observed occultations at a 
phase equal to 0.4998^ ^^. These authors re-analysed the avail- 
able photometric and spectroscopic data and finally concluded 
that <?b is likely very close to zero. 

A possible non-circular orbit might drive the dissipation of 
tidal energy. This mechanism, in turn, might be responsible for 
heating up the planet's interior and bloating its radius in con- 
sequence (Miller et al. 12009]) . The orbit of WASP-12 b is ex- 
pected to be circularised on short timescales if the tidal dissipa- 
tion constant is not significantly larger than a value typical for 
giant planets (Hebb et al. 120091 Ibgui et al. 1201 01) . Its eccentric- 
ity, if real, could be sustained by gravitational perturbations from 
an additional planet in the system (Hebb et al. 2009; Li et al. 
120101) . Such a perturbing body could affect the orbital motion of 
WASP-12 b, causing its transits to exhibit a departure from a lin- 
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Table 1. Observatories and instruments (with abbreviations identifying datasets) taking part in the campaign. 



# 


Observatory 


Telescope (abbreviation)/Instrument 


®(m) 


Mr 


1 


Observatorio del Roque de los Muchachos, La Palma (Spain) 


Nordic Optical Telescope (NOT)/ALFOSC 


2.6 


5 


2 


Calar Alto Astronomical Observatory (Spain) 


2.2-m reflector (Calar Alto)/CAFOS 


2.2 


15 


3 


National Astronomical Observatory Rozhen (Bulgaria) 


Ritchey-Chretien-Coude (2.0Rozhen)/CCD 


2.0 


10 






Cassegrain (0.6Rozhen)/CCD 


0.6 


5 


4 


Bohyunsan Optical Astronomy Observatory (South Korea) 


Bohyunsan Optical Telescope (BOT)/CCD 


1.8 


1 


5 


Osservatorio Astronomico di Bologna, Loiano (Italy) 


Cassini Telescope (Loiano)/BFOSC 


1.5 


2 


6 


Gunma Astronomical Observatory (Japan) 


1.5-m Ritchey-Chretien (Gunma)/CCD 


1.5 


2 


7 


Devasthal Observatory (India) 


Devasthal Fast Optical Telescope (DFOT)/CCD 


1.3 


1 


8 


Michael Adrian Observatory, Trebur (Germany) 


TIT (Trebur)/CCD 


1.2 


6 


9 


Mt. Lemmon Optical Astronomy Observatory (South Korea/USA) 


1.0-m reflector (Mt. Lemmon)/CCD 


1.0 


3 


10 


Konkoly Observatory (Hungary) 


Ritchey-Chretien-Coude (Konkoly)/CCD 


1.0 


2 


11 


Lulin Observatory 


Lulin One-metre Telescope (LOT)/CCD 


1.0 


1 


12 


University Observatory Jena (Germany) 


Schmidt Teleskop Kamera (Jena)/CCD 


0.9 


5 


13 


Stara Lesna Observatory (Slovakia) 


0.5-m reflector (Stara Lesna)/CCD 


0.5 


3 



Notes. Q) is the diameter of a telescope's main mirror, and N tI is the total number of complete or partial transit light curves. 



ear ephemeris (Miralda-Escude 2002; Schneider 2004; Holman 
& Murray 2005; Agol et al. l2005l Steffen et al. l2007l . 

The transit timing variation (TTV) method has been already 
used to confirm the planetary nature of systems with multi- 
ple transiting planets (e.g. Holman et al. 120101 1 and to discover 
an additional planetary companion which cannot be detected 
with other techniques (Ballard et al. 1201 j} . All these discover- 
ies are based on data acquired with the Kepler space telescope 
(Borucki et al. 1201 Oi l and, to our knowledge, no independently 
confirmed TTV signal has been detected from the ground so far] 
Maciejewski et al. d201 11 1 reported on a sign of a TTV signal 
for WASP-12 b and put upper constraints on a mass of a pos- 
sible perturber. Being encouraged by this result, we organised a 
multi-site campaign for follow-up timing observations of transits 
of WASP- 12 b. 

2. Observations and data reduction 

During three consecutive observing seasons, spanning from the 
autumn of 2009 to the winter of 2012, we acquired 61 com- 
plete or partial transit light curves. We used fourteen 0.5-2.6- 
m telescopes located in thirteen observatories distributed around 
the world. A portion of data was obtained in collaboration with 
the Young Exoplanet Transit Initiative (YETI, Neuhauser et 
al. 1201 IK A list of participating observatories and instruments 
(sorted by telescope mirror diameter) is presented in Table [T] To 
quantify the quality of each light curve, we used the photometric 
noise rate (pnr) defined by Fulton et al. (1201 II) as 



pnr : 



Vf 



(1) 



where the root mean square of the residuals, rms, is calculated 
from the light curve and a fitted model, and F is the median 
number of exposures per minute. Observations were recorded 
in Coordinated Universal Time (UTC). To exclude systematic 
errors in the recorded times, local computer clocks were verified 
with the Network Time Protocol software, usually accurate to 
better than 0.1 s. Observations were performed using the filters 
in which individual instruments were the most efficient (usually 
R). Most of telescopes were significantly defocused to reduce 
flat-fielding noise and decrease the amount of time lost due to 



CCD readout. Exposure times were kept fixed during each ob- 
serving sequence. Table [2]lists individual observations acquired 
with 2-m class telescopes while data obtained with smaller in- 
struments are given in Table [3] Short descriptions of individual 
instruments, observation details, and data reduction are given in 
Sects. 2. If 2.13 In further analysis, three amateur light curves 
(Sect. |2.14 1 were also included. 



1 Here we mean TTV detection confirmed e.g. by RV or transit sig- 
nals of the additional planet or by additional TTV observations of an 
independent group. 



2. 1 . Nordic Optical Telescope 

Five light curves were acquired in the 2011/12 observing season 
with the Nordic Optical Telescope (NOT) at the Observatorio 
del Roque de los Muchachos, La Palma (Spain), through the 
P44-102 observing programme (OPTICON 2011B/003). The 
Andalucia Faint Object Spectrograph and Camera (ALFOSC) 
was used in imaging mode, giving a field of view (FoV) of 
6'A x 6'A. The CCD was binned 2 x 2 to obtain faster readout. 
During the first run (2011 Nov 22-24) we acquired two com- 
plete light curves and one partial light curve, for which some of 
the data were lost due to bad weather. The second run (2012 Feb 
1-2) allowed us to observe two almost complete transits. The 
beginning of the event on Feb 2 was lost due to bad weather. 
Twenty-minute gaps in light curves on 2011 Nov 22 and 2012 
Feb 1 are due to the target passing through zenith, causing rapid 
motion of the field derotator. Autoguiding guaranteed that stel- 
lar centroids did not change their locations on the detector matrix 
during each night. 

The photometric measurements were performed with a 
custom-made pipeline employing differential aperture photome- 
try and standard CCD reduction procedures (de-biasing and flat- 
fielding with sky flats). Only the nearest comparison star was 
used to minimise possible flat-field systematics induced by the 
field rotator (e.g. Moehler et al. 2010). A number of apertures 
were used and the light curves with the smallest photometric 
scatter were taken as the final ones. 



2.2. Calar Alto 2.2-m telescope 

We observed 15 transits with the 2.2-m telescope at the Calar 
Alto Observatory (Spain). Most of these observations cover 
complete events and only four are partial due to bad weather. In 
the 2010/11 observing season, ten nights were awarded within 
the programme Fl 1-2.2-008. We gathered usable data in six 
nights; the other four nights were completely lost due to bad 
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Fig. 1. Light curves for transits of WASP-12 b observed with 2-m class telescopes. The residuals are shown in bottom plots. 



weather. In the 201 1/12 season, we got seven nights (programme 
F12-2. 2-009) and in all of them acquired complete transit light 
curves. Two additional light curves were obtained on 2011 Oct 
7 and Nov 10 as backup observations for the programme Hll- 
2.2-011. 

The Calar Alto Faint Object Spectrograph (CAFOS) was 
used in imaging mode. The original FoV was windowed to 
5'.6 x 5!8 and 2x2 binning was applied to shorten the read- 
out time. The telescope was autoguided during all nights except 
2012 Feb 24, when the observers were unable to use the autogu- 



ider due to significant telescope defocusing. The data reduction 



procedure was similar to that in Sect. 2.1 



2.3. 2.0-m and 0.6-m telescopes at the National Observatory 
Rozhen 

The 2.0-m Ritchey-Chretien-Coude (RCC) telescope at the 
Bulgarian National Astronomical Observatory Rozhen was used 
to observe six complete and four partial transits. During the 
2009/10 and 2010/11 seasons, a VersArray BOOB CCD cam- 
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Fig. 2. The same as Fig.^butfor 1.5-m and smaller telescopes. 
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era (1340 x 1300 pixels, 20 jim pixel size, resolution 07258 per 
pixel) with a FoV of 5!8 x 5?6 was used as a detector. The miss- 
ing part in the beginning of the event on 2009 Nov 12 as well as 
gaps on 2010 Nov 29 were caused by bad weather. 

During the 2011/12 season, the same CCD camera was at- 
tached to the telescope through the two-channel focal reducer 
FoReRo-2, giving a resolution of a.' 74 per pixel and a FoV of 
about 15' in diameter. Gaps in data on 2011 Nov 22 and 2012 
Jan 19 were caused by bad weather conditions. The beginnings 
of the events on 201 1 Dec 27, 2012 Jan 18, and 2012 Feb 23 were 
lost due to technical problems with the telescope. Red noise is 
noticeable in some of the light curves. This effect is principally 
caused by observing without autoguiding, as the telescope's au- 
toguider is unusable in the focal reducer mode. 

In addition, we acquired five transit light curves with the 
0.6-m Cassegrain photometric telescope equipped with a FLI 
PL09000 CCD camera (3056 x 3056 pixels, 12 jjm pixel size, 
resolution 07334 per pixel, FoV: 17' x 17'). 

Standard IDL procedures (adapted from DAOPHOT) were 
used for data reduction (de-biasing and flat-fielding with sky 
flats) and performing differential aperture photometry. Using the 
method of Everett & Howell (2001), several stars (4 to 6) with 
photometric precision better than 4 mmag (6 mmag for the 0.6- 
m telescope) were selected to create an artificial comparison star 
used for the differential photometry. 

2.4. Bohyunsan Optical Telescope 

The transit light curve on 2012 Dec 24 was obtained with the 
Bohyunsan Optical Telescope (BOT) at the Bohyunsan Optical 
Astronomy Observatory (BOAO). The telescope is an alt-az 1.8- 
m (f/80) reflector operated by the Korea Astronomy and Space 
Science Institute (KASI), equipped with the BOAO 4k CCD 
(15x15 /mi) imaging instrument. An available FoV is 14^6x14^6 
with a scale of 0721 per pixel. We used the 2x2 binning mode 
and windowed the FoV to increase the sampling rate. The auto- 
guider system was used and pointing errors were found to be of 
the order of a few tenths of a pixel per exposure. Observations 
were reduced in the same way as described in Sect. 2. 1 



2.6. 1.5-m Gunma Astronomical Observatory telescope 

Transits on 2009 Dec 26 and 2010 Jan 18 were observed with the 
1.5-m Ritchey-Chretien telescope at the Gunma Astronomical 
Observatory (Japan). An Andor DW432 CCD camera (1250 x 
1 152 pixels) was used as a detector providing a 12'.5 x 1 1 '.5 FoV. 

The photometric data were reduced following standard pro- 
cedures including subtraction of a median dark frame and divid- 
ing by a sky flat field. The relative magnitudes were derived with 
the aperture photometry method implemented in the customised 
software pipeline developed for the Semi-Automatic Variability 
Search sky survey (Niedzielski et al. 2003). An artificial compar- 
ison star was generated iteratively using 20-30% of the unsatu- 
rated stars with the lowest scatter. A number of aperture radii 
were tested and the aperture which produced light curves with 
the smallest scatter was finally chosen. 



2. 7. Devasthal Fast Optical Telescope 

The transit on 2011 Feb 5 was observed with the 1.3-m 
Devasthal Fast Optical Telescope (DFOT, Sagar et al. 120121 1 lo- 
cated at the Devasthal Observatory (India). A set of 2770 expo- 
sures was recorded with an exposure time of 5 s. The data re- 
duction was done as described in detail in Joshi et al. ((2009). 
To improve the signal-to-noise ratio, subsets of 20 single frames 
were co-added to produce the final light curve. 



2.8. 1.2-m Michael Adrian Observatory telescope 

The 1.2-m TIT telescope at the Michael Adrian Observatory in 
Trebur (Germany) was used to observe four and two complete 
transits in the 2010/11 and 2011/12 seasons, respectively. The 
telescope has a Cassegrain optical system and is equipped with 
a 3072 x 2048-pixel SBIG STL-6303 CCD camera (FoV: 10!0 x 
6f 7). The light curve acquired on 201 1 Jan 28 was affected by 
high airmass at the end of the run. Observations on 2011 Feb 
7 were interrupted by occasional passing clouds. Observations 



were reduced in the same way as described in Sect. 2.6 



2.5. Cassini Telescope 

A complete transit of WASP-12 b was observed on 2009 Nov 
11 and a partial one on 2011 Oct 7 with the 1.52-m Cassini 
Telescope at the Astronomical Observatory of Bologna in 
Loiano (Italy). This telescope was already successfully utilised 
to follow up several transiting planets (Southworth et al. 2010 
I2012al I2012bl Mancini et al. 120121 . It is equipped with the 
Bologna Faint Object Spectrograph & Camera (BFOSC), an in- 
strument built to allow, with a simple configuration change, the 
acquisition of both images and spectra. The detector is a back- 
illuminated EEV LN/1300-EB/1 CCD with 1300 x 1340 pixels, 
giving a FoV of 13'0 x 12f6 at a scale of 0758 per pixel. 

The observations were analysed following standard meth- 
ods. We created master bias and flat-field images by median- 
combining sets of bias images and sky flats, and used them to 
correct the science images. Aperture photometry was performed 
using the IDL/ASTROLIB implementation of DAOPHOT 
(Southworth et al. 2009). The apertures on the target star were 
placed by eye and their sizes were set to a wide range of val- 
ues to find those which gave photometry with the lowest scatter 
compared to a fitted model. Differential photometry was finally 
performed using an optimal ensemble of comparison stars. 



2.9. Mount Lemmon Optical Telescope 

Two complete (201 1 Dec 8 and 2012 Jan 1) and one partial (2012 
Jan 13) transit light curves were obtained using the 1.0-m reflec- 
tor operated by KASI at the Mount Lemmon Optical Astronomy 
Observatory (LOAO) in Arizona, USA. For all observations we 
used the ARC 4k CCD instrument in the 2x2 binning mode. The 
CCD has a pixel size of 15 x 15 /im, a scale of 0741 per pixel, 
and a FoV of 28' x 28'. During the transits the weather was clear 
but not photometric with partial clouds early in the evenings and 
later clearing up during the runs. Observations were reduced in 



the same way as described in Sect. 2.6 



2.10. Konkoly Observatory 

Transits on 201 1 Jan 27 and 201 1 Dec 5 were observed with the 
1 .0-m RCC Telescope at the Piszkesteto Station of the Konkoly 
Observatory (Hungary), equipped with a 1340 x 1300-pixel PI 
VersArray 1300b NTE CCD (FoV: 7' x 7', image scale a/32 per 
pixel). Flat-fielding using sky flats, dark correction, and aperture 
photometry was performed with standard IRAF routines and a 
dedicated gnu-r script as described in Szabo et al. ( 120 10b . 
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2.11. Lulin One-meter Telescope 

A transit on 2010 Feb 10 was observed with the Lulin One- 
meter Telescope (LOT) at the Lulin Observatory, operated by 
the National Central University of Taiwan. A VersArray:1300B 
CCD camera with a FoV of 11' x 11' was used as a detec- 
tor. Observations were reduced in the same way as described 
in Sect. ESI 

2. 12. University Observatory Jena 

Five light curves were acquired with the 0.9/0. 6-m Schmidt tele- 
scope at the University Observatory Jena in GroBschwabhausen 
near Jena (Germany). The CCD-imager Schmidt Teleskop 
Kamera (STK, Mugrauer & Berthold [20T0l > has a 52'.$ X 52'.$ 
FoV which allows numerous comparison stars to be observed 
simultaneously. 

Data acquired on 2011 Feb 8, 2011 Mar 2 and 2011 Mar 3 
were reduced as described in Sect. |2.6| The remaining the nights 
(201 1 Nov 21 and 22) were reduced with the IRAF routines after 
subtracting a dark frames and applying flat-field correction with 
sky flats. After applying aperture photometry, differential mag- 
nitudes were calculated using the algorithm described in Broeg 
et al. (12005b . This algorithm calculates an artificial standard star 
out of all stars in the field, weighted by their magnitude errors, 
and compares the brightness of each star, including our target 
star WASP-12, against this standard star. 

2.13. Stara Lesna Observatory 

Three light curves were obtained at the Stara Lesna observa- 
tory in Slovakia. Observations were performed with the 0.5-m 
Newtonian telescope, on which an SBIG ST10 MXE CCD cam- 
era with 2184x1472 6.8-yum pixels was mounted. The scale was 
0756 per pixel, corresponding to a FoV of 24' x 16'. 

The light curve obtained on 2010 Mar 9 exhibits correlated 
noise because observations were acquired without autoguiding. 
Standard correction procedures (bias, dark, and flat-fielding) and 
then aperture photometry was performed with the C-munipack 
packag^] To generate an artificial comparison star, at least 20- 
30% of stars with the lowest photometric scatter were selected 
iteratively from the field stars brighter than 2.5-3 mag below 
saturation level. To measure instrumental magnitudes, various 
aperture radii were used. The aperture which gave the smallest 
magnitude scatter was applied to generate a final light curve. 

2. 14. Amateur light curves 

Light curves for transits on 2011 Oct 30 and 2012 Jan 7 were 
acquired at Castelnuovo Magra (Italy) with a Meade LX200 
GPS 10 inch telescope equipped with an f/6.3 focal reducer and 
an SBIG ST8 XME CCD camera. Science frames were taken 
through the Baader Yellow 495 Longpass filter with an exposure 
time of 300 s. The CCD frames were acquired with the MaxIm 
DL 4 software and then reduced with MaxIm DL 5. The data re- 
duction procedure included removing dark frames, flat-fielding, 
and aperture photometry against nearby bright stars. 

Observations of the transit on 2011 Dec 27 were acquired 
with a CCD camera based on a Sony ICX429ALL sensor, 
mounted on an 0.2-m Newtonian telescope located near Lucca 
(Italy). The data were also collected through the Baader Yellow 
495 passband filter and an exposure time of 260 s was used. The 

2 http://c-munipack.sourceforge.net 
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Fig. 3. Additional amateur light curves obtained with 20-cm 
reflecting telescopes. The mid-transit times were found to be 
BJDtdb 2455865.6130± 0.0012 (2011 Oct 30), 2455923.4609 ± 
0.0026 (2011 Dec 27), and 2455934.3725 ± 0.0009 (2012 Jan 
07). 

equatorial fork mount of the telescope was guided by a modi- 
fied Webcam coupled with a 60 mm telephoto lens. Images were 
acquired by the AstroArt v4 software and a standard data re- 
duction procedure including aperture photometry was performed 
with the MaxIm DL 5 programme. To generate an artificial com- 
parison star, eight reference stars of magnitudes close to that of 
WASP-12 were used. 

2.15. Light curve pre-processing 

Even the best-quality light curves may be affected by photomet- 
ric trends caused, for instance, by colour difference between the 
target and comparison stars and changes of airmass during a run. 
If they are not removed properly they may seriously affect transit 
timing. To remove trends, a model that consists of a transit sig- 
nal adopted from Maciejewski et al. d201U and a second-order 
polynomial was fitted to each light curve analysed in this work. 
We used the jktebop code (Southworth et al. I2004al I2004bt . 
which is based on the ebop program (Etzel fTWTt Popper & Etzel 
I1981I ) and allows photometric trends to be approximated with a 
polynomial up to order five. A best-fitting trend was subtracted 
from each light curve and then magnitudes were transformed 
into fluxes normalised to unit out-of-transit flux level. 

3. Results 

3. 1 . System parameters 

The Transit Analysis Package^] (tap, version 2.104, Gazak et al. 
120121 ) was used to analyse transit light curves and refine sys- 
tem parameters. TAP uses the MCMC (Markov Chain Monte 
Carlo) approach, with the Metropolis-Hastings algorithm and a 
Gibbs sampler, to determine model parameters. Wavelet-based 
techniques are employed to account for correlated noise while 
estimating error values. It has been shown that this approach 
provides the most reliable parameter results and error estimates 
(e.g. Hoyer et al. 120121 1. The code employs the quadratic limb- 
darkening (LD) law to model the distribution of the flux on a 
stellar disc. 

To refine the system parameters, the best-quality light curves 
acquired with 2-m class telescopes were selected for a simul- 
taneous fit. The selection criteria included sub-millimagnitude 
precision of a dataset (pnr < 1 .0 mmag) and transit phase com- 
pleteness higher than 75%, including at least an ingress or/and 
egress preceded or followed by out-of-transit observations and 
a flat-bottom phase. Nineteen light curves were found to ful- 
fil these criteria, including observations on 2010 Feb 02 from 

3 http://ifa.hawaii.edu/users/zgazak/IfA/TAP.html 
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Table 4. Parameters of the WASP-12 system from the joint anal- 
ysis of 19 sub-millimagnitude-precision light curves and the re- 
fined linear ephemeris. 



Value 



Parameter 

Orbital inclination, i b 
Scalled semimajor-axis, ab/R, 
Planetary to stellar radii ratio, Rb/R, 
Transit parameter, b 
Linear LD coefficient, u 
Quadratic LD coefficient, v 
Observed planetary radius R b (Rj up ) 

Effective planetary radius, Rf (R lup ) 1 .90 ± 0.09 

Stellar density, p, (p Q ) 0.315 ± 0.007 

Orbital period, P b (d) 1.0914209 ± 0.0000002 

Cycle-zero transit time, T Q (BJD TD b) 2454508.97718 ± 0.00022 



82°96+°o 50 

o.ujo_ 0021 
H 73+ 00005 
355 +bm 

301 +0 (H| 
n ?4 +0.<K 
"•^ -0.08 

1.86 ±0.09 



Maciejewski et al. (201 1 ). We ran ten MCMC chains, each con- 
taining 10 6 steps. The individual chains were combined to get 
final posteriori probability distributions. The following param- 
eters were linked together for all light curves: orbital period 
Pb, orbital inclination semimajor-axis scaled by stellar ra- 
dius ab/R*, planetary to stellar radii ratio Rb/R*, and linear (m) 
and quadratic (v) LD coefficients. The mid-transit times, airmass 
slopes, and flux offsets were allowed to vary separately for in- 
dividual light curves. A circular orbit was assumed. An analy- 
sis of relations between the parameters reveals significant cor- 
relation or anti-correlations (with the Pearson correlation coeffi- 
cients ranging from 0.71 to 0.96) between it, and ab/R*, i\, and 
Rb/R*, ab/R* and Rb/R*, Rb/R* and v, and u and v. The best- 
fitting parameters were determined by taking the median value 
of marginalised posteriori probability distributions which were 
found to be unimodal. The 15.9, and 85.1 percentile values of 
these distributions were taken as upper and lower 1 cr errors. 
The results are presented in Table |4] where the impact param- 
eter (defined as b — |r cos ib) is also given. The final model of 
WASP-12 b's transit signature is plotted in Fig. [4] 

The LD coefficients of the best-fitting model were compared 
to the theoretical values calculated from the tables from Claret 
& Bloemen d20TTT>, and derived with the EXOFAST appleQ 
(Eastman et al. 1201 21) . The coefficients were bilinearly interpo- 
lated from the tables assuming the host star effective temperature 
and surface gravity as given by Hebb et al. (J2009 ). The theo- 
retical coefficients, u t = 0.288 and v t = 0.317 agrees with the 
derived values (u = 0.301+°^i and v _ 0.24+J}^) within 1 cr. 

To check for any evolution of the system parameters, for 
instance as a result of a periastron rotation, the fitting proce- 
dure was repeated with ib, ab/R*, and Rb/R* allowed to vary 
between individual light curves. The linear and quadratic LD 
coefficients were allowed to vary around values derived in the 
previous run, under the Gaussian penalty defined by the derived 
errors. We found that none of these parameters exhibited a statis- 
tically significant variation. In addition, we compared individual 
mid-transit times derived in both approaches, and no discrepancy 
greater than 1 cr was found. 

The planet WASP-12 b is known to depart significantly from 
a spherical shape. The ratio of its polar and sub-stellar-point radii 
was found to be ~ 1.14 (Budai l201 11 1. The observed planetary ra- 
dius Rb is the cross section of the planet at inferior conjunction. 
It can be calculated from ab/R* and Rb/R* assuming a value of 
the semi-major axis of Ob from Hebb et al. (2009). The effective 
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Fig. 4. Best-fitting model of WASP-12 b's transit (red line) laid 
over the 19 best-quality, phase-folded light curves. 



planetary radius, Rf, which is more representative of the physi- 
cal parameters of the planet, is greater. Following the numerical 
results of Budaj (EoTTT) . we get Rf = 1.90 ± 0.09 fl Jup . This 
value is identical to 1 .90 ± 0.09 Rj up determined by Maciejewski 
et al. (120111) and within 1-cr of the results of Chan et al. (1201 11 1 
and Southworth (12012I I. The mean planetary density was found 
to be pb = 0.206 + 0.043 pj up . This value is smaller than found 
by previous studies, because we have accounted for the effec- 
tive planetary radius in our calculations, but is consistent with 
them within 1 cr. The surface gravitational acceleration we find, 
gb = 9.8 ± 0.4 m s~ 2 , is also noticeably smaller than previous 
studies, for the same reason. 

A faint stellar companion to WASP-12 has been recently 
discovered at a separation of ~1" (Bergfors et al. 12012b . It is 
~4 mag fainter in the i' band and its spectral type is in the region 
of K0 to M0. As it is not resolved in any of our observations, its 
light contributes to the total flux of the system. This makes the 
transit depth shallower by 0.2-0.3 mmag, the precise amount 
depending on its spectral type. Due to this, we estimate the final 
planetary radius to be greater by 0.6-1.0%, hence pb and gb will 
be slightly smaller than the values given above. 

The stellar parameters agree with previous studies. We find 
the stellar radius to be R* = 1 .63 + 0.07 R Q (assuming the value 
of Ob from Hebb et al. 120091 1. its mean density to be p* = 0.315 ± 
0.007 p Q , and its surface gravity to equal logg* = 4.14 ± 0.03 
(cgs units). 

3.2. Transit timing 

The individual transit times were determined by fitting a tem- 
plate transit light curve with tap. The fitting procedure started 
with initial parameters whose values are taken from the best- 
Table [4j. They were allowed to vary 



fitting model (Sect. 3.1 



http://astroutils.astronomy.ohio-state.edu/exofast/limbdark.shtml 



in the MCMC analysis under Gaussian penalties whose scales 
were determined by the parameters' errors. This approach al- 
lowed best-model uncertainties to be included in a mid-transit 
time budget of each light curve. Ten chains of length of 10 5 
steps were used for each light curve. The first 10% of itera- 
tions of each chain were discarded before calculating the pa- 
rameter values and errors. The mid-transit times, converted 
to Barycentric Julian Dates in Barycentric Dynamical Time 
(BJDtdb, see Eastman et al. 120101 1. and their errors are given 
in Tables [2] and [3] 

In addition, we re-analysed two datasets from Maciejewski 
et al. (2011) and obtained mid-transit times 2455230.40669 ± 
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Fig. 5. O-C diagram for WASP-12 b's transit timing. The filled circles and squares denote our observations and literature ones, 
respectively, with errors smaller than 40 s. The open symbols - circles and squares for our and literature data, respectively - mark 
points with timing errors greater than 40 s. The open triangles denote the amateur times included in this work. The postulated 
sinusoidal variation is sketched with a red line. Individual observing seasons, covered by our high-precision observations, are 
zoomed-in in the bottom. The literature data are taken from Hebb et al. (2009), Chan et al. (2011), Maciejewski et al. (1201 11) . Sada 
et al. d20T2l , and Cowan et al. d2012] ). 



0.00019 and 2455254.41871 ± 0.00045 BJD TD b- These values 
are well within error bars compared to the original mid-transit 
times of 2455230.40673 ±0.000 1 1 and 2455254.4 1 887+0.000 1 4 
BJDtdb, respectively. We also re-analysed a high-accuracy light 
curve acquired with the 2-m Liverpool Telescope and published 
in Hebb et al. ( 120091 ). The mid-transit time for this early epoch 
(no. 6) was found to be 2454515.52496 ± 0.00043 BJD TD b as- 
suming times in original data are in HJDutcJJI The original mid- 
transit time was not published in Hebb et al. (2009 ) but it is listed 
in Campo et al. d20TT1) as 2454515.52542 ± 0.00016 BJDtdb- 
Both determinations differ by less than 1 cr. 

Mid-transit times of a single transit observed with different 
telescopes provide an empirical test of the timing accuracy and 
error estimates. We collected two or more light curves for twelve 
transits. Eight of them were observed with two telescopes and 
four with three telescopes. The median difference between indi- 
vidual times was found to be 0.4 cr and never exceeded 1 cr. This 
finding supports the reliability and consistency of our timing sur- 
vey. Furthermore, we compared mid-transit times derived with 
tap and jktebop. The latter code offers Monte Carlo (MC) simu- 
lations, a bootstrapping algorithm, and the residual-shift (prayer- 
bead) method to estimate the uncertainties of parameters. The 



5 It is stated that the light curve is in ID^-rc (L. Hebb 2012, priv. 
comm.). However, the mid-transit time derived by us and the value 
published in Campo et al. ( 201 1 1 clearly indicate that the original light 
curve was converted into HIDutc- 



largest of the errors from the methods were conservatively taken. 
Both codes give consistent results with the median difference be- 
tween the determined times to midpoint equal to 13 s or 0.13 cr. 
The median of the ratios of the tap to the jktebop timing errors 
was found to be 1.33. This value includes partial light curves for 
which the tap code usually gives errors greater by a factor of up 
to four. 

A linear function was fitted to the transit times, to de- 
termine the orbital period and the time of transit at cy- 
cle zero, To. Cycle zero was taken to be the epoch given by 
Hebb et al. ([2009). We fitted the transit times from new and 
literature data by the method of least-squares. The data point 
for epoch zero was skipped because it comes from a global 
fit based on a number of individual light curves (Hebb et al. 
120091 We found T = 2454508.97718 + 0.00022 BJD TD b and 
P h = 1 .09 14209 ± 0.0000002 d. The individual mid-transit errors 
were used as weights during the fitting. The reduced ;f 2 is 1.29, 
which corresponds to a p-value of 0.06. Hence a null hypothesis 
assuming constant P\, could not be rejected. The O-C (Observed 
minus Calculated) diagram for transit timing is plotted in Fig. [5] 

More interestingly, the same conclusion cannot be drawn if 
the mid-transit time sample is limited to the best-quality data. 
We now consider only those mid-transit times which have uncer- 
tainties below 40 s (i.e. a half of the TTV amplitude suggested 
in Maciejewski et al. 1201 II) . The linear fit was repeated, result- 
ing inx 2 = 1-99 and a /rvalue of 1.6 x 10 -3 . This result allows 
the null hypothesis to be questioned. Moreover, a coherent pat- 
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Fig. 6. Periodogram of the O-C residuals (continuous black line) 
and spectral window (dashed red line, multiplied by 10 for clar- 
ity). The dashed vertical line shows the empirical FAP level of 
1%. The arrow marks the most significant peak which, as can be 
clearly seen, is not caused by the structure of gaps in the time 
domain. The shortest period investigated is given by the Nyquist 
frequency. 



tern in the O-C diagram is noticeable (Fig. [5]). In the next step, a 
model consisting of a linear trend and a sinusoidal variation was 
fitted. After removing the linear component, a Lomb-Scargle pe- 
riodogram (Lomb 1976 Scargle 1982 ) reveals a peak at a period 
of ~510 epochs (Fig. |6j. The false alarm probability (FAP) of 
this signal is 0.7% and was determined empirically by a boot- 
strap resampling method. The procedure randomly permutes the 
O-C values retaining the original observing epochs and then cal- 
culates the Lomb-Scargle periodogram to determine the power 
of the strongest signal. After 10 5 such trials, the FAP of the sig- 
nal in the original dataset is taken to be the fraction of resam- 
pled periodograms that produced higher power than the original 
dataset. From the least-squares fit, we find that the possible TTV 
signal has a period of 500 + 20 epochs and a semi-amplitude of 
0.00068 ± 0.00013 d (Fig. [5]). The errors are taken from the co- 
variance matrix of the fit. The reduced^- 2 of this model is 0.91. 



3.3. Two-planet model 

A TTV signal shown by the transiting planet may reflect gravi- 
tational perturbations driven by an undetected additional planet 
in the system (e.g. Ballard et al. 2011). To check what plane- 
tary configurations may produce timing variations with the pos- 
tulated periodicity, three-body simulations were performed with 
the mercury package (Chambers 1999) employing the Bulirsch- 
Stoer integrator. A hypothetical additional planet of a probe mass 
of 40 M @ was put on orbits with semimajor axes ranging from 
0.007 to 0.068 AU (orbital periods between 0.22 and 5.5 d, re- 
spectively) with a step equal to 10~ 6 AU. The semimajor axis 
of WASP-12 b and the mass of the central star were set to the 
values found by Hebb et al. (2009). The eccentricity of the per- 
turber, e c , was varied from 0.0 to 0.3 in steps of 0.05. Scenarios 
with circular and eccentric orbits (e\, equal to 0.02, 0.03 and 
0.05) of the transiting planet were considered. The simulations 
assumed coplanar and edge-on orbits. The calculations covered 
1500 orbital periods of the transiting planet, i.e. slightly greater 
than the time span of the timing observations. Lomb-Scargle pe- 
riodograms were calculated for the synthetic O-C diagrams, to 
determine the dominant periodicity. Numerous solutions that re- 
produce the postulated TTV signal were identified. 

To alleviate the solution degeneracy, the publicly available 
RV measurements were reanalysed with the Systemic Console 
software (Meschiari et al. 12009b . We used data from Hebb et al. 



(2009), Husnoo et al. d20TTb , and Albrecht et al. (120121) . The 
datapoints in transit phase (+1.6 h from predicted mid-transit 
times) were removed to avoid interference from the Rossiter- 
McLaughlin (RM) effect (Rossiter [l924l McLaughlin \WM that 
manifests itself as a distortion of the spectral lines of the host 
star during transit. Data published by Hebb et al. ([2009) were 
collected in the 2007/8 observing season with the SOPHIE spec- 
trograph coupled with the 1 .9-m telescope at the Observatoire de 
Haute Provence. We selected 17 RV measurements with errors 
between 9 and 13 m s for further analysis. Husnoo et al. (1201 11 1 
presented further data acquired with SOPHIE in the 2008/9 and 
2009/10 seasons. The spectra acquired on 2009 Mar 21 and 2009 
Apr 1 3 have half the signal-to-noise ratio and noticeably greater 
errors than the other observations. This could indicate that these 
two exposures were affected by cloud, which, in turn, could in- 
troduce a systematic shift in the RV measurement (N. Husnoo 
2012, priv. comm.). We noticed that one of these lower-quality 
points is an outlier at the 3 cr level regardless of the RV model un- 
der consideration. These two points were therefore not included 
in any further analysis. It is worth noting that keeping these out- 
lying points does not quantitatively change the ensuing results 
but noticeably degrades the quality of the fits to the RVs. We also 
rejected measurements from 2009 Jan 17 when the signatures of 
the RM effect were observed. We finally used four and six dat- 
apoints from the 2008/9 and 2009/10 seasons, respectively. RVs 
reported by Albrecht et al. (120121 1 were acquired in the 201 1/12 
season with the HIRES spectrograph on the Keck telescope. We 
used 13 measurements in our study. 

The RV dataset from Husnoo et al. (201 1) was split into two 
separate subsets according to observing season. We noticed that 
allowing a variation in the RV offsets between individual observ- 
ing seasons results in significantly better fits. This effect may be 
caused by the low number of datapoints, the influence of addi- 
tional bodies on long-term orbits in the system, or instrumental 
effects. The latter is a justifiable possibility because SOPHIE, 
when operated in the High Efficiency mode, has been found to 
produce deviations up to several dozens of m s from a velocity 
zero point on a timescale of a few months (Husnoo et al. 1201 11 1. 

In our analysis, the RV signature of tides raised in the host 
star by the transiting planet is not considered because its ex- 
pected amplitude of 4.8 m s _1 (Arras et al. I20121 > is smaller than 
the typical errors of the RV measurements. 

The periodogram of the RV residuals versus the best-fitting 
model for WASP-12 b on a circular Keplerian orbit (Fig. |7]i 
shows a signal at approximately 3.6 d. Its empirical FAP was 
found to be 0.08%, from 10 5 trials with the bootstrap resampling 
method. A preliminary Newtonian RV modelling indicate that 
the residuals may be caused by an additional ~0. 1 Mj up planet 
on a non-circular orbit with an eccentricity <? c of between 0.2 
and 0.3. Intriguingly, our three-body simulations performed with 
the Mercury code show that a system configuration for a 0.1 
Mj up perturber may reproduce the postulated timing variations 
of WASP-12 b if e c is close to 0.3 and is between 0.03 and 
0.05. 

To model the RV and transit timing datasets simultaneously, 
we used the Systemic Console with the Bulirsch-Stoer integra- 
tor and a precision parameter of 10~ 14 . The times of RV observa- 
tions were transformed to BJDtdb- To put additional constraints 
on the orbital parameters of WASP-12 b, the transit dataset was 
enhanced with five mid-occultation times from Campo et al. 
(1201 11 1 and Croll et al. (1201 II) . The light-travel time correction of 
22.8 s was taken into account. The Lopez-Morales (2010) mid- 
occultation point was not considered in the final iterations be- 
cause of its much greater uncertainty and possible systematic 
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Fig. 7. Periodogram of RV residuals for WASP-12 b on a 
Keplerian orbit (continuous black line) and the spectral window 
(dashed red line, multiplied by 10 for clarity). The arrow indi- 
cates the strongest peak at 3.58 d that is not associated with sig- 
nificant power in the spectral window. The shortest period is set 
by the radius of the star. The dashed and dotted horizontal lines 
show empirical FAP levels of 1% and 0.1%, respectively. 
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timing offset (see Campo et al. l2011l and Croll et al. 201 1 for dis- 
cussion). The orbital periods, masses, eccentricities, longitudes 
of periastron, and mean anomalies were allowed to vary for both 
planets. The offsets for individual RV subsets were also free pa- 
rameters. The best-fitting model was found iteratively with the 
simulated annealing and Nelder-Mead minimisation algorithms. 
The joint x 2 statistic accounts for the contributions from the RV 
and transit time datasets, x\v an d^ 2 r respectively: 
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x 2 = 



1 



N RV + N ti - N 6t 



[(Mw - iVfitklv + A X 1] 



(2) 



where Mw, Mr, and A% t are the numbers of RV measurements, 
mid-transit times, and fitted parameters, respectively, and A is an 
arbitrary weight (Meschiari et al. 2009). The x\y of the best- 
fitting model is approximately five times greater than x 2 r , which 
suggests that the RV errors are underestimated. To balance the 
significance of both data sets, A = 5 was specified in the final 
iterations. The MCMC method was used to determine parameter 
uncertainties. The MCMC chain was 10 6 steps long, and the first 
10 s configurations were discarded. The scale parameters were 
set empirically in a series of attempts to get the acceptance rate 
of the MCMC procedure close to the optimal value of 0.25. For 
each parameter the standard deviation was taken as the final error 
estimate. The best-fitting parameters of the two-planet model are 
given in Table[5] Fig.[8]shows the RV components from WASP- 
12 b and the postulated planet, as well as the modelled TTV 
signal. 

The RV datasets from the 2008/9 and 2009/10 seasons were 
found to be shifted against the 2007/8 dataset by -16.4 + 3.3 
and +1.1 + 2.9 m s _1 , respectively. The value of the joint x 2 was 
found to be 3.0 with individual components xh, and;^ equal to 
2.5 and 0.4, respectively. For comparison, a single-planet model 
with a circular and eccentric orbit have joint x 2 values of 7.9 
and 6.3, respectively. In the two-planet model, the final rms of 
the RV dataset was reduced to 7.4 m s _1 from 12.5 m s _1 for 
a single -planet model with WASP-12 b on a circular orbit, and 
from 9.6 m s _1 for a single-planet model with WASP-12 b on an 
eccentric orbit. 



3.4. Dynamical stability 

The relatively large eccentricity of the postulated planet WASP- 
12 c raises the question of whether the system is dynamically 




600 800 
Epoch 



1000 1200 1400 



Fig. 8. Joint RV and transit timing model for the WASP-12 sys- 
tem including two planets. Upper, the RV variation produced by 
WASP-12 b. Middle: the RV signal generated by the postulated 
planet WASP-12 c. Bottom: the O-C diagram for the transit tim- 
ing of WASP-12 b with a modelled TTV signal and the residuals 
shifted by -0.0026 d for clarity. 



stable. In this paper, we consider only Newtonian, point-mass 
mutual interactions as a first-order model of the dynamics of the 
system. The planet WASP-12 b has extremely short orbital pe- 
riod, so it is not sufficient to study the long-term evolution of the 
system with close-in companion WASP-12 c based only on the 
conservative dynamics. Due to the close proximity of WASP- 
12 b to its host star, a full dynamical model would need to ac- 
count for general relativistic effects, rotational and tidal defor- 
mations, and the tidal dissipation of energy (e.g. Mardling & 
Lin [2002; Migaszewski & Gozdziewski 120091 . If these effects 
are considered, the orbital geometry and long-term evolution of 
the whole, mutually interacting system may be strongly affected 
since its formation, possible migration, and tidal evolution. We 
postpone this complex problem to another work. 

A cursory inspection of the putative system parameters 
(Table [5} reveals that it is close to the 10b:3c mean-motion 



10 



G. Maciejewski et al.: Multi-site campaign for transit timing variations of WASP-12 b 
Table 5. Parameters of the two-planet WASP-12 system from the joint analysis of RVs and transit timing. 



Parameter Planet b Planet c 



Orbital period, P (d) 


1.0914281 ± 


: 0.0000015 


3.58484 ± 


0.00011 


Minimum mass, in sin ( (Mj up ) 


1.356 ± 


: 0.009 


0.098 ± 


0.009 


Mean anomaly" , A(°) 


310.61 


±0.25 


184.6: 


± 1.1 


Eccentricity, e 


0.0447 ± 


: 0.0043 


0.284 ± 


0.002 


Longitude of periastron, a) (°) 


274.44 


±0.03 


223.8 : 


±2.5 


Semimajor axis, a (AU) 


0.0229 ± 


: 0.0008 


0.0507 ± 


0.0018 


Mass* , m (M Jup ) 


1.366 + 


: 0.011 


0.098 ± 


0.009 



{a) given for the initial epoch BJD TDB 2454509.3871 



assuming a coplanar system with the inclination derived from transit light curve modelling. 



resonance (MMR^j This justifies our interest in the short- 
term dynamics of the system. We reconstruct the structure of 
the phase-space in terms of the maximal Lyapunov exponent 
(Benettin et al. 1 19761 1 which is expressed through the so called 
fast indicator Mean Exponential Growth factor of Nearby Orbits 
(MEGNO, see e.g Cinco tta & Simo 120001 Cincotta et al. 120031 
Gozdziewski et al. 2001 for details). This dynamical character- 
istic helps us to classify given sets of initial conditions as regular 
(leading to quasi-periodic evolution of the system, stable over in- 
finite period of time) or chaotic (leading to irregular phase-space 
trajectories). The later might lead to geometric changes of the 
orbits rendering the system unstable in the short-term timescale 
of 10 4 -10 5 outer periods (P c ). By the design, the MEGNO indi- 
cator when converged to 2 over a given integration time, guar- 
antees orders of magnitude longer time-scale of the dynamical 
(Lagrange) stability. However, one must be aware that we con- 
sider a relatively narrow time window of the system evolution, 
restricted to the conservative, Newtonian point-masses dynam- 
ics. 

To visualise the orbital stability and dynamical neighbour- 
hood of the WASP-12 system, we computed a high-resolution 
dynamical map in the (a c ,e c )-plane, see Fig.[9} with the help of 
our new CPU cluster software mechanic (Slonina et al. I2012I I. 
Each point (initial condition) was integrated for 10 s P c . Such a 
timescale is long enough to detect even weak, high-order MMRs. 
Two MMRs (56c: 17b and 33c: 10b), which are the closest match 
to the nominal position of the system (marked with the asterisk 
symbol), are labelled on the scan. Thanks to the high resolution 
of the dynamical map, other, even weaker MMRs are also iden- 
tifiable. The 79c:24b MMR and the 89b:27c MMR (both not la- 
belled), which are between the 23b:7c MMR and 56c: 17b MMR, 
are the closest to the nominal position of the system. Bearing in 
mind very small error range of the semimajor axis of the outer 
planet in the dynamical map, ~ 10~ 6 AU, (see an explanation 
below), we found that it might encompass these resonances. Yet 
they are very narrow and weak due to their high order. There is 
only a small probability that the system might be trapped in such 
resonances. The error range of eccentricity safely separates the 
WASP-12 system from the zone of global chaos, visible roughly 
above e c ~ 0.32. This zone appears due to proximity to the or- 
bits' crossing zone and overlapping of the MMRs. 

Note that the dynamical map in Fig. [9] is computed for a 
particular, nominal mass of the parent star and all orbital ele- 
ments, besides those ones constituting the map coordinates, are 
fixed at their best fit values. In general, for a fixed nominal mass 
of the star, the semi-major axis derived as osculating element 
at some initial epoch is essentially influenced by uncertainties 



of the planetary masses and orbital periods in accord with the 
Hid Kepler law: 



k\M t + m c )P, 
(4^) 2 



2\(l/3) 



(3) 



as the most uncertain parameters (here k is the Gaussian gravita- 
tional constant), though the planetary and stellar mass relation is 
indirect due to the iV-body observational model. Therefore, the 
formal error of a c in the dynamical map should be considered as 
weakly dependent on the error of the (fixed) stellar mass. 

If the stellar mass is varied, the relative positions of the mean 
motion resonances do not change, as they depend essentially on 
the mass ratios 



— = 1 + — 

M, + m h \ M„ 



I 1 + 



nib 



(4) 



6 The Xc:Yb notation means the ratio of orbital periods of the inner 
and the outer planet, respectively. 



and the ratios of the orbital periods. Hence, the dynamical map 
might be considered to be representative for different masses of 
the parent star. We also emphasize, following Michtchenko & 
Malhotra (2004), that in the limit of small planetary masses and 
eccentricities, the Newtonian, conservative three-body dynamics 
depend on the mass ratios and semi-major axes, and not on the 
individual masses and semi-major axes as the system seems to 
be well separated from low-order MMRs (see also Migaszewski 
& Gozdziewski 120091 . 



4. Discussion 

The two-planet model predicts that the eccentricity of WASP- 
12 b might be non-zero. This, in turn, might affect the system 
parameters. The fitting procedure was repeated with non-zero e^ 
and a>b treated as fixed parameters, values of which were taken 
from Table [5] The noticeable difference was found for z'j, and 
at>/R* which were found to be degrees and 2.93 ±0.02, 

respectively. Compared to the circular-orbit scenario, the ;'b and 
flb/7?* determinations are smaller by 0.5 <x and 2.5 <x, respec- 
tively. Adopting these values in the calculations results in 
and R* greater by 4% and 2%, respectively, hence smaller values 
for pb, gb, p*, and logg*. As the non-zero eccentricity of WASP- 
12 b is speculative, we settle for the scenario with the circular 
orbit in Sect.lXTl 

Statistical studies of hot Jupiter candidates in the Kepler 
sample show that these planets are rare in multi-transiting sys- 
tems (Latham et al. 1201 11 1. The occurrence rate is below 5%. 
The solitariness of hot Jupiters is postulated to be a result of 
the dynamical evolution of their planetary systems. On the other 
hand, there are systems with hot Jupiters accompanied by non- 
transiting planets on wide orbits, e.g. HAT-P-13 (Bakos et al. 
120091 . HAT-P-17 (Howard et al. l20T2l . and Qatar-2 (Bryan et 



11 



G. Maciejewski et al.: Multi-site campaign for transit timing variations of WASP-12 b 



<Y> index 




0.0505 0.0507 0.0509 

semi-major axis a c [au] 



0.0511 



Fig. 9. Dynamical map of the WASP-12 system in the (a c , em- 
plane, computed in terms of the MEGNO fast indicator. These 
two orbital parameters are varied across the map, while other or- 
bital elements and masses are kept at their nominal values from 
Table[5] The star symbol, located at a c = 0.050681 AU, indicates 
the nominal system. Colours encode the dynamical stability of 
the system: yellow is for strongly unstable configurations, and 
blue is for stable, quasi-periodic systems. A few relevant MMRs 
are labelled. The raw resolution of this map is 1440 x 900 data- 
points. Each initial condition at the map has been integrated for 
10 5 outer orbits. 



al. 12012b . These additional planets have been discovered with 
the RV technique. Their gravitational influence on the transiting 
planets is too small to be detected through TTVs. Interestingly, 
single hot or warm Jupiter candidates which exhibit TTV sig- 
nals have been found in the Kepler sample (Szabo et al. I2Q121 >. 
The frequency of hot Jupiter candidates with periodic TTVs is 
estimated to be 18%. More than half of these positive detections 
may be of the perturbative nature that allows scenarios with an 
additional close planet or exomoon to be considered. 

The ground-based timing of transiting hot Jupiters has 
brought no positive results so far, though preliminary detections 
have been announced. The characteristics of a spurious TTV sig- 
nal for HAT-P-13 b are the most similar to the case of WASP- 
12 b. Pal et al. d20 lit observed a 0.01 d departure from a lin- 
ear ephemeris. This finding was confirmed by Nascimbeni et 
al. d20 lit who proposed a sinusoidal variation with a period of 
~1 150 d and semi-amplitude of 0.005 d. Such a signal could not 
be caused by the known HAT-P-13 c planet, suggesting the exis- 
tence of a third planetary component. Then, Fulton et al. d201 11 1 
showed that all mid-transit times except the one from an ear- 
lier work by Szabo et al. (2010) are consistent with a linear 
ephemeris. A similar conclusion was reached by Southworth et 
al. (2012c) who found no evidence for the existence of periodic 
timing variations. 

Bearing in mind the case of HAT-P-13 b, we observed var- 
ious phases of the TTV signal postulated for WASP-12 b: a 
maximum (the 2009/10 season), two minima (the 2010/11 and 
2011/12 seasons), and portions of ascent and descent trends. 
Furthermore, we selected mid-transit times of the highest qual- 
ity for the final analysis. They come from different instruments 
and data reduction pipelines, thus systematic errors are expected 
to be minimised. The timing scatter between consecutive transits 
(up to four in the 201 1/12 season) and runs (separated from each 
other by at least 10-12 epochs) was found to be much smaller 
than the amplitude of the postulated long-term variation. 

Near-ultraviolet (near-UV) transmission spectroscopy dur- 
ing transits of WASP-12 b shows that the planet has an exo- 



sphere extending beyond the planet's Roche lobe (Fossati et al. 
1201 Obi ). Moreover, these data show that transits begin earlier and 
are longer than at optical wavelengths. This phenomenon may 
be interpreted as a fingerprint of an accretion stream escaping 
from the planet through the inner Lagrangian point (Lai et al. 
2010). It is also explicable as the effect of a bow shock induced 
by the orbital motion of the planet with a magnetic field in the 
stellar coronal material (Vidotto et al. 2010; Llama et al. 1201 II ). 
A characteristic of such dynamical interactions is temporal vari- 
ation implied by the variable magnetic field of the central star 
(Vidotto et al. 120 1 lb . The changes in phases of the early ingress 
would affect the transit duration and hence the time of a mid- 
transit. This, in turn, would produce a TTV. The effect of the 
early ingress, however, is not expected to be detected in the op- 
tical bands. In Sect. 3.1 we analysed only complete transits, and 
found no variation in the transit duration. We also compared the 
values of the rms in ingress, flat bottom, and egress phases of 
these light curves after phase folding. Changes in the transit du- 
ration would smear the ingress and egress data points in phase, 
causing a greater rms. We found the rms values of 0.80, 0.76, 
and 0.77 mmag for ingress, flat bottom, and egress, respectively. 
These values are compatible with each other, so again we find 
no sign of transit duration variations (TDV). 

Other mechanisms which generate TTVs may also be re- 
jected. An effect induced by an exomoon is expected to produce 
a number of harmonic frequencies in the periodogram of the 
timing residuals (Kipping 2009). The interactions with a distant 
stellar or planetary component on a wide orbit would generate a 
long-period TTV through the light-travel-time effect (Montalto 
2010). The amplitude of the TTV caused by the Applegate mech- 
anism (Applegate 1992), which invokes quasi-periodic varia- 
tions in the quadrupole moment of magnetically active stars, 
is expected to be a few seconds over the course of ~10 years 
(Watson & Marsh 2010). The effects of transit parallax (Scharf 
120071) and stellar proper motion (Rafikov 2009 ) are to small to 
be detected for the WASP-12 system at present. 

Croll et al. (12011b note that e\, is poorly constrained if a>b is 
close to 90° or 270°, and this is the case for WASP-12 b. The 
literature determinations of this angle range from 270° (Croll et 
al. l20TT1 ) to 286° (Hebb et al. 120091 . Our two-planet model gives 
0Jb = 274° 44 ± 0°03. We obtained values of e b cos u h = 0.0035 ± 
0.0002 and e h sin w b = -0.0446 + 0.0034, and both of them are 
clearly far from zero. For the latter parameter, the consistency 
with the literature determinations is between 1.1 and 1.3 cr for 
Campo et al. ( 120111 and Croll et al. d2011b . respectively. In the 
case of eb cos u^, the values are seemingly discrepant if they are 
compared to our value obtained for the initial epoch 2454509 
BJD. 

Our model neglects non-Keplerian sources of the apsidal 
precession including general relativistic effects, the quadrupole 
field caused by rotational flattening of the star and planet, and 
tidal potentials of the star and planet. This approach is justi- 
fied because these effects are usually small, and hence generate 
precession rates much lower than the precession due to the sec- 
ond planet. However, it has been shown that the precession rate 
due to the tidal potential of WASP-12 b may be non-negligible 
(Ragozzine & Wolf 2009 ). In this case, the precession period is 
predicted to be several years or decades, depending on the re- 
sponse of the planetary body to a tidal potential. 

WASP-12 is a quiet star with no signatures of a structured 
magnetic field (Fossati et al. 120 1 Oat and a chromospheric ac- 
tivity index log R' HK = -5.5 (Knutson et al. 120101 Krejcova & 
Budaj 12012b . To estimate the expected amount of stellar RV jit- 
ter generated by the non-uniform convection (Saar et al. 1998), 
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we analysed the results of Isaacson & Fischer (2010) who stud- 
ied RV stability for more than 2600 main sequence and sub- 
giant stars with spectral types between F5 and M4. We selected 
a sample of 310 stars whose colour index B - V and activity in- 
dex \ogR' HK are close to these of WASP-12. As WASP-12 has 
B — V — 0.6 mag we considered stars with 0.5 < B — V < 0.7 
mag, and with log R' HK < -5.0. The median jitter was found 
to be 2.6 ms" 1 . The stellar jitter of WASP-12 is too high in a 
single-planet model, being 12.5 and 6.7 m s in the cases of cir- 
cular and eccentric orbits for WASP-12 b, respectively. In the 
two-planet model, the jitter was found to be 2.7 ms _1 , which is 
consistent with the expected value from the Isaacson & Fischer 
(120101 1 survey. To verify the jitter of WASP-12, we used the 
Albrecht et al. (2012) RV dataset which contains precise mea- 
surements from a single night. In such a short time interval of a 
few hours, the RV variation is generated by stellar jitter, planets, 
and possible instrumental trends. The latter two components can 
be temporarily approximated with a single Keplerian orbit. This 
approach yields a jitter of 2.7 ms~' for WASP-12, in agreement 
with the expected value. 

The radii of known transiting planets with masses close to 
0.1 Mj up range from 0.21 Rj up (Kepler-22 b, Borucki et al. 12012b 
to 0.73 flj U p (Kepler-35(AB) b, Welsh et al . 120121 . If the postu- 
lated planet transits the host star, the corresponding flux drop is 
between 0.02% and 0.2%. If the planet is bloated to 1 Rj up , the 
expected transit depth is 0.4%, so such a signal could be detected 
with a 1-m class telescope. The geometric transit probability for 
an eccentric orbit (Kane & von Braun 2008 ) of the postulated 
planeQwas found to be 13 ± 4%, i.e. almost three times lower 
than for WASP-12 b. If the system is coplanar, the additional 
planet could produce grazing transits. 

The co-planarity of the system was found to be supported 
by numerical experiment in which the mutual inclination of both 
planets, A;, was varied from to 30° with a step size of 2° for 
Ai < 10° and 5° for 10° < Ai < 30°. The goodness of the fit 
for the two-planet model was found to decrease with increasing 
Ai. Scenarios with Ai < 10° result in < 1, so the two-planet 
WASP-12 system is likely to be coplanar. This finding testifies 
in favour of quiescent mechanisms by which planets are trans- 
ported inwards with the planets' orbital orientation preserved 
(i.e. interactions between planets and protoplanetary discs). On 
the other hand, the obliquity of the host star (the sky-projected 
angle between the stellar spin and a planetary orbital axis) was 
found to be 59+H degrees (Albrecht et al. 120121 ). Such a high 
value suggests that WASP-12 b migrated inwards by a process 
changing the relative orientation of the stellar and planetary axes 
(i.e. planet-planet interactions). It cannot be completely ruled 
out that the axial misalignment has a primordial origin unrelated 
to the evolution of the planetary system. It has been shown that 
an interaction between a magnetosphere of a young magnetic 
star and circumstellar disc may produce a significant spin-orbit 
misalignment (Lai et al. |201 11 ). 



5. Conclusions 

Transit timings of WASP-12 b and precise velocities of its host 
star give marginally significant support for the presence of the 
second planet, when both datasets are considered separately. 
However, a joint analysis results in a two-planet model which 
better explains observations than single-planet scenarios. This 
approach removes the degeneracy of solutions which is typical 



7 Here we assume a planetary radius of 0.7 Rj up . 



of the TTV method. Our finding allows us to advance the hy- 
pothesis that WASP-12 b is not the only planet in the WASP-12 
system, and that there is an additional 0.1 Mj up body on a 3.6-d 
period eccentric orbit. The dynamical analysis of the two-planet 
system shows that it is located in a relatively wide stable zone. 
This reinforces our hypothesis. However, we note that the puta- 
tive planetary companion remains under the threshold of reliable 
detection. Thus, further photometric and spectroscopic follow- 
up observations of high precision are needed to shed new light 
on the architecture of the WASP-12 system. 
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Table 2. Transit light curves observed with the >1.8-m telescopes. 



# 


Date UT 


Epoch 


Telescope 


Filter 


r 


pnr 
(mmag) 


Tmid (BJD TDB ) 

2450000+ 


O-C 

(d) 


HQ 


1 


2009 Nov 12 


586 


2.0Rozhen 


-^Cousins 


4.60 


0.78 


5148.5514™ 


+0.0016 


— 


2 


2010 Nov 29 


936 


2.0Rozhen 


^Cousins 


1.77 


1.40 


5530.54705™? 
554? 55271 +uuuu28 

-0.00029 


-0.00009 


— 


3 


2010 Dec 11 


947 


2.0Rozhen 


-^Cousins 


1.77 


0.79 


-0.00007 


yes 


4 


2011 Jan 03 


968 


Calar Alto 


-^Cousins 


0.83 


0.70 


5565.4718+°$;| 


-0.0008 


— 


5 


2011 Jan 04 


969 


Calar Alto 


^Cousins 


1.05 


0.71 


5566.56333!° S 

i- y nn nnn t 4-M (Will >(-, 

5600 39791 +U UUU2b 
5601.48963™ 
5612 4042 

JU,i,wi -0.00 16 


-0.00070 


yes 


6 


2011 Feb 07 


1000 


Calar Alto 


^Cousins 


0.98 


0.76 


-0.00018 


yes 


7 


2011 Feb 08 


1001 


Calar Alto 


^Cousins 


0.83 


1.06 


+0.00013 


— 


8 


2011 Feb 19 


1011 


Calar Alto 


^Cousins 


0.94 


1.11 


+0.0005 


— 


9 


2011 Mar 02 


1021 


Calar Alto 


-^Cousins 


0.88 


0.75 


5623.3178!°™ 
5842 69343+ 000091 


-0.0001 




10 


2011 Oct 07 


1222 


Calar Alto 


^Cousins 


0.82 


1.05 


-0.00010 


— 


1 1 


2011 Nov 10 


1253 


Calar Alto 


-^Cousins 


0.84 


0.99 


5876.52839! 00 
5887.44183+ ; 

COQQ CTTT 1 +0.00037 

3888.^3 /1_ . 00038 


+0.00081 


yes 


12 


2011 Nov 21 


1263 


2.0Rozhen 


^Cousins 


1.54 


0.86 


+0.00004 


yes 


13 


2011 Nov 22 


1264 


NOT 


^Bessel 


1.72 


0.85 


+0.00050 


yes 


14 


2011 Nov 22 


1264 


2.0Rozhen 


-^Cousins 


1.54 


1.27 


i- i-i n n i-^i c i J_(l 11/111 w 

000049 


+0.00030 


— 


15 


2011 Nov 23 


1265 


NOT 


^Bessel 


1.82 


0.57 


5889.6253+°; 00 ' 


+0.0006 


— 


16 


2011 Nov 24 


1266 


NOT 


^Bessel 


1.20 


0.75 


5890.71638!;* 


+0.00033 


yes 


17 


2011 Dec 24 


1293 


BOT 


^Cousins 


0.44 


1.74 


5920 18500 + °™ 5 

5923 45846 
j;/ij.tjotu 0M)20 


+0.00059 


— 


18 


2011 Dec 27 


1296 


2.0Rozhen 


^Cousins 


1.23 


0.62 


-0.00021 


yes 


19 


2012 Jan 18 


1316 


2.0Rozhen 


^Cousins 


1.23 


0.82 


I- r-\ a r- i~i i-i s i-r-i J-( \ IW 11170 

5945 28657 


-0.00053 


yes 


20 


2012 Jan 19 


1317 


Calar Alto 


^Cousins 


0.83 


1.16 


5946 37835+ °036 
5946 37S72 +UUUUJ4 

J7tU.J/0/Z, _n.pOQ33 


-0.00016 


— 


21 


2012 Jan 19 


1317 


2.0Rozhen 


-^Cousins 


1.23 


1.00 


+0.00021 


— 


22 


2012 Jan 20 


1318 


Calar Alto 


-^Cousins 


0.82 


0.87 


5947 46973+ ™2i 


-0.00020 


yes 


23 


2012 Jan 21 


1319 


Calar Alto 


^Cousins 


0.82 


1.10 


5948 56075+° 00028 

— Q QQQ28 

cnrn 00018 

5959.47548 „„„„ „ 
5960 56687+° ™° 


-0.00060 


— 


24 


2012 tea 01 


1329 


NOT 


^Bessel 


1.54 


0.51 


-0.00009 


yes 


25 


2012 Feb 02 


1330 


NOT 


^Bessel 


1.58 


0.45 


-0.00012 


yes 


26 


2012 Feb 12 


1339 


Calar Alto 


-^Cousins 


0.83 


0.72 


5970.38980!°;° 00 | 


+0.00003 


yes 


27 


2012 Feb 13 


1340 


Calar Alto 


-^Cousins 


0.83 


0.89 


5971.48128!°;° oo p 


+0.00009 


yes 


28 


2012 Feb 23 


1349 


2.0Rozhen 


^Cousins 


1.23 


0.85 


5981 30398+° ™ 4 


-0.00000 


yes 


29 


2012 Feb 24 


1350 


Calar Alto 


-^Cousins 


0.88 


0.87 


SQ87 QQS7?+°™047 

5983.48682!; 
6005.31541+™ 


-0.00018 


yes 


30 


2012 Feb 25 


1351 


Calar Alto 


^Cousins 


0.87 


0.94 


+0.00000 


yes 


31 


2012 Mar 18 


1371 


2.0Rozhen 


^Cousins 


1.54 


0.87 


+0.00017 


yes 



Notes. T is the median number of exposures per minute, pnr is the photometric noise rate, r m y is the mid-transit time, O-C is the difference from 
the linear ephemeris given in Sect. [3/2] and HQ flags light curves used in the joint analysis (Sect. [3~T) . 
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Table 3. Transit light curves observed with the <1.8-m telescopes. 



# 


Date UT 


Epoch 


Telescope 


Filter 


r 


pnr 

(mmag) 


T m id (BJD TDB ) 
2450000+ 


O-C 

(d) 


1 


2009 Nov 1 1 


585 


Loiano 


Gunn-r 


0.49 


1.04 


5147.45861™ 


+0.00020 


2 


2009 Nov 22 


595 


0.6Rozhen 


^Bessel 


0.49 


2.91 


5158 ' 3734 -oE, 


+0.0008 


3 


2009 Nov 23 


596 


0.6Rozhen 


^Bessel 


0.49 


2.89 


5159 46377+™ 95 

ju;».tuj/ '_n 00088 


-0.00026 


4 


2009 Nov 24 


597 


0.6Rozhen 


^Bessel 


0.49 


3.10 


5160 55549+™™ 


+0.00003 


5 


2009 Dec 26 


626 


Gunma 


-^Cousins 


1.00 


1.93 


5192 20540+™-*' 


-0.00127 


6 


2010 Jan 18 


647 


Gunma 


^Cousins 


1.00 


1.81 


5215.12869™ 


+0.00218 


7 


2010 Feb 10 


668 


LOT 


^Bessel 


1.15 


2.22 


5238.04586+° : 000 g 


-0.00048 


8 


2010 Mar 09 


693 


Stara Lesna 


^Cousins 


0.97 


2.58 


5265 3313+ 0015 


-0.0006 


9 


2010 Oct 24 


903 


Trebur 


^Bcsscl 


1.23 


1.70 


5494 530 18+™ 48 

Jt^t.jjuio _ 00047 


-0.00007 


10 


2010 Oct 24 


903 


0.6Rozhen 


^Besscl 


0.65 


2.37 


5494 53061+° 00063 


+0.00035 


11 


2011 Jan 13 


977 


0.6Rozhen 


^Bcsscl 


0.65 


2.99 


5575 29703+ 00075 

JJ/J.z.7/ —0.00079 


+0.00163 


12 


2011 Jan 27 


990 


Konkoly 


-^Cousins 


0.95 


1.63 


5589 48296+™' 1 

Q 00094 

5590 57597+o mioo 

/ J" / —0.00092 


-0.00091 


13 


2011 Jan 28 


991 


Trebur 


^Bcsscl 


0.87 


2.52 


+0.00067 


14 


2011 Feb 05 


998 


DFOT 


"Cousins 


0.57 


1.14 


5598 21556+ 00036 
5600 39720+™" 


+0.00032 


15 


2011 Feb 07 


1000 


Trebur 


^Bessel 


0.87 


1.90 


-0.00088 


16 


2011 Feb 07 


1000 


Stara Lesna 


-^Cousins 


1.51 


2.28 


5600 3978+° 0010 

JUUU.J7/0 _o oio 


-0.0003 


17 


2011 Feb 08 


1001 


Trebur 


^Bessel 


0.88 


2.01 


5601 48999+° 00086 

-0.00088 


+0.00048 


18 


2011 Feb 08 


1001 


Jena 


^Bcsscl 


1.04 


2.03 


5601.49023! ;g 


+0.00073 


19 


2011 Mar 02 


1021 


Jena 


^Bcsscl 


1.25 


2.15 


5623.31864+° 000 « 


+0.00072 


20 


2011 Mar 03 


1022 


Jena 


^Bcsscl 


0.82 


2.77 


5624 41007+° 00086 

juz-t.tiuu/ a.00082 


+0.00073 


21 


2011 Mar 03 


1022 


Stara Lesna 


^Cousins 


1.48 


2.86 


5624 4104+° 0011 

JUZ-t.tlUt n nnj . 

5842.6942+°°°; 3 


+0.0010 


22 


2011 Oct 07 


1222 


Loiano 


Gunn-r 


0.45 


1.79 


+0.0007 


23 


2011 Nov 21 


1263 


Trebur 


^Bcsscl 


0.88 


1.95 


5887.44184!° 000 « 
5887.44237.™ 


+0.00006 


24 


2011 Nov 21 


1263 


Jena 


^Besscl 


2.04 


1.83 


+0.00059 


25 


2011 Nov 22 


1264 


Jena 


^Bcsscl 


2.57 


1.90 


COQQ C1T)C + (1. 00(157 

DoS5.Djj/D_ 00005 „ 

5901 63002+™™ 

— Q.QQQ83 


+0.00005 


26 


2011 Dec 05 


1276 


Konkoly 


^Cousins 


0.33 


2.58 


-0.00023 


27 


2011 Dec 08 


1278 


Mt. Lemmon 


-^Cousins 


0.82 


2.42 


5903 81351+™' 3 

J7UJ.01JJ1 „„„„„, 

5926 73339+° ™- 

' J J J7 -QJ)0089 


+0.00041 


28 


2011 Dec 31 


1299 


Mt. Lemmon 


^Cousins 


0.82 


1.67 


+0.00046 


29 


2012 Feb 14 


1311 


Mt. Lemmon 


^Cousins 


0.82 


1.70 


5939.8286!° °°! 8 


-0.0014 


30 


2012 Mar 19 


1372 


Trebur 


^Bcsscl 


0.88 


2.03 


6006.40642!°;°™° 


-0.00025 



Notes. T is the median number of exposures per minute, pnr is the photometric noise rate, 7__<i is the mid-transit time, and O - C is the difference 
from the linear ephemeris given in Sect. |3.2| 



